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Abstract. This paper is concerned with the initial-boundary value problem 
for the Einstein equations in a first-order generalized harmonic formulation. 
We impose boundary conditions that preserve the constraints and control the 
incoming gravitational radiation by prescribing data for the incoming fields 
of the Weyl tensor. High-frequency perturbations about any given spacetime 
(including a shift vector with subluminal normal component) are analyzed using 
the Fourier-Laplace technique. We show that the system is boundary-stable. In 
addition, we develop a criterion that can be used to detect weak instabilities with 
polynomial time dependence, and we show that our system does not suffer from 
such instabilities. A numerical robust stability test supports our claim that the 
initial-boundary value problem is most likely to be well-posed even if nonzero 
initial and source data are included. 
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1. Introduction 

Most attempts to solve the Einstein equations numerically are based on the Cauchy, 
or 3 + 1 formulation of general relativity, in which one foliates spacetime into three- 
dimensional spacelike hypersurfaces of constant time t. Initial data satisfying the 
constraint equations are prescribed on the t = surface and the evolution equations 
are integrated towards the future. Typically, one truncates the domain of integration 
and only aims to find solutions on a compact spatial manifold with artificial timelike 
boundaries, on which boundary conditions must be specified. These should satisfy a 
number of requirements: 

(i) The initial-boundary value problem (IB VP) should be well-posed, i.e., for given 
initial and boundary data a unique solution should exist and it should depend 
continuously on the data. 

(ii) The boundary conditions must be compatible with the constraints in the sense 
that if the constraints are satisfied initially then they remain satisfied at all times. 

(iii) The boundary conditions should (in some sense) control the incoming 
gravitational radiation. 

The last requirement is of particular importance if one aims to obtain reliable 
information about the gravitational radiation emitted by compact sources such as 
coalescing binary black holes. One way to describe the incoming radiation is in terms 
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of the incoming fields of the evolution system that the Weyl tensor obeys by virtue of 
the Bianchi identies. These incoming fields are proportional to the Newman-Penrose 
scalar tyy, a quantity that is invariant under infinitesimal coordinate transformations 
of Kerr spacetime. If the outer boundary is placed sufficiently far away from the source 
then one expects 



to be a reasonable "no-incoming-radiation" boundary condition. Conditions of this 
type have been considered by a number of authors Q] [2] [31 Q] • 

Since the initial-boundary problem for Einstein's equations was first studied by 
Stewart [S] with numerical relativity in mind, this has been a very active field of 
research. However, there is currently only one formulation, due to Friedrich and Nagy 
PP, that satisfies all the above requirements for the fully nonlinear vacuum Einstein 
equations. It uses a tetrad formalism and evolves the components of the Weyl tensor as 
separate variables. All the boundary conditions are written in maximally dissipative 
form and well-posedness of the IBVP follows from standard theorems for symmetric 
hyperbolic systems with such boundary conditions [HI [S] ■ It would be desirable 
to obtain a similar result for a metric formulation of the Einstein equations which 
does not evolve the Weyl tensor separately, such as most of the formulations currently 
used in numerical relativity (see however [5] for an implementation of a tetrad-based 
formulation). In some cases, well-posedness has been proved (at least partially) for 
boundary conditions that satisfy conditions (i) and (ii) but not (iii) (see for example 
HOI EH E2I ESI EH EH US] , in particular harmonic formulations [UJ Ql QU EU] ) . The 
major obstruction is the fact that the physical boundary conditions (JIJ as well as 
(in general) the constraint-preserving boundary conditions contain second rather than 
first derivatives of the metric (such boundary conditions are said to be of differential 
type) . Hence they are not in maximally dissipative form and the standard theorems do 
not apply. Progress can still be made in a certain approximation (typically in the high- 
frequency limit or for linearizations about flat space) using Fourier-Laplace transform 
methods [^[55]. This technique yields rather strong necessary conditions for well- 
posedness, although sufficient conditions are not known for boundary conditions that 
arc of differential type. Applications to the Einstein equations and related problems 
can be found in [SJ E3 E3] E] ESI 120] • I n particular, boundary conditions satisfying 
both requirements (ii) and (iii) were analyzed in 0] I25| . Necessary conditions for 
well-posedness were verified but numerical experiments indicated that the systems 
still suffered from instabilities. Using rather different methods based on semigroup 
theory, Nagy and Sarbach (23 have recently proved well-posedness of the IBVP (with 
radiation-controlling boundary conditions) for the linearized Einstein equations in the 
ADM formulation in a certain gauge. The gauge condition arises from a minimization 
principle and involves a fourth-order elliptic equation for the lapse function (which 
might be expensive to solve numerically). 

In this article, we consider a symmetric hyperbolic formulation of the Einstein 
equations based on generalized harmonic (GH) coordinates (see [27] and references 
therein). Such coordinates x a obey the scalar wave equation 



with a source H a that may depend on the coordinates and the spacetime metric ij) a f, 
(but not on its derivatives). In this gauge, the Einstein equations reduce to a system 
of coupled nonlinear wave equations, with principal parts 



*o = 



(1) 



Ox a = H' 



a 



(2) 



Olp ab 



~ 0. 



(3) 
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Recently, harmonic coordinates have played an important role in the first successful 
simulations of the inspiral, merger and ringdown of binary black holes by Pretorius 
I2EH2H]- In [27]; a first-order formulation of this system was derived, with particular 
emphasis on controlling the constraints in numerical evolutions. Boundary conditions 
satisfying requirements (ii) and (iii) above were constructed, and a partial result on 
the well-posedness of these boundary conditions was stated. The purpose of this 
paper is to elaborate on and generalize this result and to provide both analytical and 
numerical evidence that the system at hand is a strong candidate for a well-posed 
initial-boundary value formulation of the Einstein equations. 

In section [2] we review the evolution equations and boundary conditions derived 
in [22] • The Fourier-Laplace analysis of the initial-boundary value problem is carried 
out in section [2] Whereas the well-posedness result in [22] assumed that the shift 
vector was tangential at the boundary, we lift that restriction here and allow for an 
arbitrary shift. Because several characteristic fields propagate along the normal lines of 
the spacetime foliation, this may cause additional characteristic fields to be incoming, 
which complicates the analysis. Nevertheless, the result stated in [22] carries over 
to the general case. (We still assume for the analysis that the normal component of 
the shift is subluminal at the boundary, i.e., we disregard situations where either all 
the modes are outgoing or all the modes are ingoing. If all the modes are outgoing, 
e.g. inside a black hole, no boundary conditions are to be imposed. See sectionQj]and 
|22] for numerical tests of the GH evolution system involving this situation.) 

The usual necessary condition for well-posedness in the Fourier-Laplace 
framework amounts to verifying that a certain complex determinant has no zeros 
with positive real part. We strengthen this result by showing that our system also 
obeys the Kreiss condition |21|. which is stronger and implies that the solution can 
be estimated in terms of the boundary data. Related concepts of well-posedness are 
briefly discussed. It has been claimed |241 0] that even if the Kreiss condition is 
satisfied, weak instabilities with milder than exponential time dependence might exist 
if nontrivial initial data are included. In section 01 we show rather generally that 
the Kreiss condition excludes any instabilities with polynomial time dependence. We 
also demonstrate that for a bad choice of gauge boundary conditions that do not 
satisfy the Kreiss condition, a weak instability does exist. In order to supplement our 
analytical results, we perform a series of numerical robust stability tests in section 
This involves adding random data to the initial and boundary conditions and to the 
right-hand-sides of the evolution equations. The background spacetime is taken to be 
either flat space (including a shift) or Schwarzschild. A summary of the results and 
conclusions are given in section [B] 

2. The initial-boundary value problem 

We begin by presenting the initial-boundary value problem for the Einstein equations 
in generalized harmonic gauge that is considered in this paper. This section follows 
closely [22] • Some additional details concerning the constraint-preserving boundary 
conditions in the case of a non-tangential shift at the boundary are provided. 

2.1. Evolution equations 

We consider the first-order formulation of the generalized harmonic evolution 
equations ^ derived in |27| . The fundamental variables are the spacetime metric 
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ip ab and its first derivatives <&i ab = diipab and n a ;, = —t c d c 4> a b- Here i a oc d a t denotes 
the unit timelike normal to the t = const, hypersurfaces. Lower-case Latin indices 
from the beginning of the alphabet denote 4-dimensional spacetime quantities, while 
lower-case Latin indices from the middle of the alphabet are spatial. 
The evolution equations take the form 

d t ^ab^(l+ll)N k d k ^ab, (4) 

d t U ab ~ N k d k U ab - Ng k *d k <S> lab + 7i72^V fc cW a6 , (5) 

dt$i ab ~ N k d k <$> lab - NO t U ab + Nj 2 d^ ab , (6) 

where ~ indicates that only the principal parts of the equations are displayed. The 
spatial metric gij , lapse function N and shift vector N l are related to the 4- metric via 
the standard ADM form of the line element, 

ds 2 = ^ ab dx a dx b = -N 2 dt 2 + g i3 \dx l + N l dt)(dx j + N j dt). (7) 

From now on we will choose 71 = —1, which ensures that the evolution equations 
are linearly degenerate. The parameter 72 was introduced in |27j in order to damp 
violations of the constraint 

C iab = ditpab - $iab = 0. (8) 

There is also a parameter 70 hidden in the source terms, which is designed to damp 
violations of the generalized harmonic gauge constraint J2Jl, based on a suggestion in 
|3(J| . It will play no role in this discussion. 

The system is symmetric hyperbolic for any choice of parameters 71 and 72 . 
Hence the Cauchy problem is well-posed 22 . The characteristic variables and speeds 
in the direction of a unit spacelike vector ni are given by 

u° ab = ipab, speed 0, (9) 

ulf = U ab ± $„ afe - 72^6, speed -N n ±N, (10) 
u Aab = ®Aab, speed — N n , (11) 

where here and in the following, an index n denotes contraction with n, (e.g., 
v n = UiV 1 ) and upper-case Latin indices denote projection with = g^ — n^ij 
(e.g., v A = PaiV 1 ). 

We consider a finite spatial domain of integration £1 with smooth boundary dCl. 
The spatial coordinate location of the boundary remains fixed in time, i.e., d/dt is 
tangential to the 3-dimensional timelike boundary dQ x [0, 00). Boundary conditions 
have to be prescribed for the incoming characteristic fields, where now m refers to 
the outward-pointing unit spacelike normal to dCl. Note that the number of incoming 
fields depends on the value of the normal component N n of the shift at the boundary. 
For —N < N n < 0, only the fields u]^ are incoming; for < N n < N, the fields 
and u\ ab are incoming. 

2.2. Constraint-preserving boundary conditions 

The boundary conditions should ensure that if the constraints vanish initially then 
they vanish at all times. In other words, the subsidiary evolution system that the 
constraints obey as a consequence of the main evolution equations has to be well- 
posed, with the unique solution being the trivial one. 

In our case, the primary constraints are the harmonic gauge constraint @, which 
can be written in the form 

C a = V b T abc + H a = (12) 
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O^abc being the Christoffel symbols of the metric lpab)^ an d the definition constraint 
Ciab ©• If we define the first-order constraint variables 

Cijab = 29[i$j] a b = 2d[iCj] ab , (13) 

T a = -t c d c C a + . . . , (14) 

Cia = d,C a - g jk C ljka + IgJ^Cijcd + ■■■ (15) 

(omitting terms proportional to Ci ab ), then the constraint evolution system can be 
written in the simple form |27] 

dtC a =s 0, (16) 

d t Ta ^ N%F a + Ng ij diC ja , (17) 

d t C ia ~N'd 3 C ia + Nd l T a , (18) 

d t C iab ~ 0, (19) 

d t C ijab ~ N k d k C ljab . (20) 

This system is clearly symmetric hyperbolic, and the characteristic variables and 
speeds are 



r 0± — T znC 


speed - TV" ± TV, 


(21) 


c 1 -C 


speed 0, 


(22) 


c Aa = CAa, 


speed -iV", 


(23) 


c iab = Ciafc, 


speed 0, 


(24) 


c ya6 — ^ijab. 


speed -iV". 


(25) 



Consider first the case —N< N n ^ 0. The only incoming constraint fields are 
c a~ ■ We impose completely absorbing boundary conditions on them: 

eg" = (26) 

(= denoting equality at the boundary). These can be translated into conditions on 
the normal derivatives of four of the main incoming fields uzZ by noting that 

c °- ~ V2 [k^ d \ - \k a ^ cd \ d n ulj 

+±p A a yj cd d A n cd - p AB d A n Ba - p AB t b d A $ Bba (27) 

+\t a $ cd P AB d A $Bcd + l2(P AB d A yj Ba - \P A a r d d A ^cd) 

+ P AB d A $ nB a ~ \P A a4> Cd d A <$> ncd , 

where we have introduced the ingoing null vector k a = (t a - n a )/V2. 

In other words, the constraint-preserving boundary conditions imply boundary 
conditions on the normal derivatives of the following projection of u ab , 

\P ab P cd - 2l (a P b) < c fc d ) + lj b k c k d ] u\- , (28) 



p(C)cd 1- 
r ab a cd 



where we have also introduced the outgoing null vector l a = (t a + n a )/^/2. This 
projection has four degrees of freedom. Of the remaining six degrees of freedom of 
u ab ' t w0 will be fixed by the physical boundary conditions and four by the gauge 
boundary conditions. 

Next, we consider the case < N n ^ N. Now the fields c 2 Aa and cfj ab are 
incoming as well. In addition to (|26|l . we impose the boundary conditions 

= CnAbc - dnU Abc - d A <&nbc, (29) 
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which are conditions on the normal derivatives of the main characteristic fields 
u \ab- The remaining incoming constraint fields c Aa and c\ Bab also need boundary 
conditions. However, together with the physical and gauge boundary conditions, 
we have already used up all the incoming modes of the main evolution system and 
cannot impose any further boundary conditions actively. Fortunately, it turns out 
that c Aa = and c ABab = as a consequence of the constraint-preserving boundary 
conditions we have already imposed and the evolution equations. We prove this 
important point in the appendix, using the Fourier-Laplace method (we only consider 
the limit of high-frequency perturbations about a fixed background spacetime in that 
part of the analysis). 

To summarize, we have imposed homogeneous maximally dissipative boundary 
conditions for the constraint evolution system. Because of general theorems for 
quasilinear symmetric hyperbolic systems with such boundary conditions 13 IEJ- 
it follows that the constraint evolution system is well-posed and the unique solution 
is the trivial one. 

Setting to zero the incoming fields of the constraint evolution system at the 
boundary as in (j26|) is the standard prescription used in most works on constraint- 
preserving boundary conditions (e.g., |S1 EH Q3] EH E])- This ensures that any 
constraint violations generated in a numerical evolution leave the computational 
domain without reflections (at least for normal incidence). In harmonic formulations, 
one sometimes considers the simpler alternative C a = |20| , which does not involve any 
derivatives of the fundamental fields. These Dirichlet boundary conditions are clearly 
also consistent with the constraints but they constitute a reflective boundary condition 
for the constraint evolution system (fToU2THl . Numerical tests of such boundary 
conditions for an axisymmetric version |31| of the Z4 system |32| (which is very similar 
|3"U] to the generalized harmonic evolution system considered here) indicate that they 
can cause significant reflections of constraint violations |25|. 

2.3. Physical boundary conditions 

The physical boundary conditions used in [27J are designed to control the gravitational 
radiation entering the domain through the boundary, following the prescription used in 
mUOUE]- As explained in the introduction, the incoming radiation can be described 
in terms of the ingoing characteristic fields w~ b of the Weyl tensor evolution system, 

w~ b = 2(P a c P b d - \P ab P cd )k e kfC ced} . (30) 

(The fields w~ b are proportional to the Newman-Penrose scalar \&o for a null tetrad 
containing the null vectors k a and l a defined above.) As a boundary condition, we 
impose 

< b = d t h%. (31) 

(p) 

The data h ab may be used to inject a gravitational wave through the boundary, 
for instance. Note that the expression for w~ b is only unique up to multiples of the 
constraints Cij a b related to the index ordering in the first-order reduction. For a 
particular choice of index ordering, we can write w~ b in the form 

ur b ~ (P a A P b B - \P ab P AB ) [d n {u l ^ B + l2 u° AB ) 

+ LV-7M-,,/; - //'"VM-//.., - P CD dc$DAB (32) 

+t C d A <£nBc ~ d A IlBn - t C d A $>Bnc] ■ 
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Hence (|31H provides a boundary condition on the normal derivatives of 

€ )cd <d = (Pa c P b d lP ab P cd Ka- (33) 

In |33| , a hierarchy of boundary conditions has recently been developed that are 
perfectly absorbing for linearized gravitational radiation up to some arbitrary angular 
momentum number L. They can be viewed as successive improvements of the ^0 = 
boundary conditions and take the form 

d t (r 2 l a d a ) L - l {r^ Q )\ r=R = 0, (34) 

where r is an areal radial coordinate and it is assumed that the boundary is a sphere 
of radius R. We remark here that the stability result derived in this paper carries over 
to these improved boundary conditions (cf. section T3.2|l . 



2.4- Gauge boundary conditions 

The remaining components of u7 correspond to gauge degrees of freedom. To 
understand this, one should note that the generalized harmonic condition does not 
fix the coordinates completely: there is still a remaining gauge freedom x a — -> x a + £ a 
for infinitesimal coordinate displacements £ a that satisfy the wave equation. We may 
exploit this gauge freedom in order to choose the four remaining boundary conditions 
in any way we like. Our viewpoint here is that the gauge boundary conditions should 
guarantee that the IBVP is well-posed, while still admitting arbitrary boundary data 
in order to be able to impose a variety of gauge conditions. 

(G) 1 — 

The simplest way appears to prescribe data h ab for those components of u ab 
that are not specified by the constraint-preserving and physical boundary conditions, 



p(G)cd 1- 
r ab U cd 



xc sd p(C)cd p (P)cd 
°(a°b) - ^ab ~ ^ab 



cd 



= [2l {a k b) l c k d + k a k b l c l d - 2k [a P b) c l d ] u l c - = h§\ (35) 

Equivalcntly, we can write (|3 51) in the simpler form 

l b ul b =h^. (36) 

These (with M G > = 0) are the gauge boundary conditions used in the numerical 
implementation of \27\. 

For illustrational purposes only, we shall also consider the following alternative set 
of gauge boundary conditions, which will turn out to be ill-posed (to varying extent): 

t h U ab = hf\ (37) 

To linear order, these amount to prescribing (the time derivative of) the lapse and 
shift at the boundary. 



3. Fourier-Laplace analysis 

In this section, we discuss necessary conditions for wcll-posedness using the Fourier- 
Laplace technique. We begin with a review of the general theory, mainly following 
|SJ |2] and highlighting some open questions. We then apply it to the generalized 
harmonic evolution system for the case of high-frequency perturbations about any 
arbitrary spacetime. 



Stable radiation- controlling boundary conditions 8 
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Consider a linear symmetric hyperbolic system of evolution equations 

d t u = A l diU, (38) 

where it is an m-dimensional vector, the A 1 are m x m constant symmetric matrices, 
and 1 ^ i ^ n. We take the spatial domain of integration to be the quarter-space 
t ^ 0, x 1 ^ 0, — oo < x 2 , x 3 , . . . , x n < oo. Initial data are prescribed at t = 0, 

u(0,x i ) = f(x i ). (39) 

Let I denote the number of negative eigenvalues of A 1 (i.e., the number of incoming 
modes). We impose boundary conditions at x 1 = of the form 

S%u = d t h, (40) 

where the S l are constant I x m matrices and ft, is a vector of boundary data (this 
non-standard form of the boundary conditions is required for our application in section 

We allow for the boundary to be (uniformly) characteristic, which is the case for 
most applications in physics, including the system discussed in this paper. After a 
suitable orthogonal transformation of the variables, we may write 

(«) 

where is the k x k zero matrix and A is an (m — k) x (m — k) matrix. Similarly, 
we split 

u = (z,v) T . (42) 

The initial-boundary value problem l|38H40|l can be solved by means of a Laplace 
transform with respect to time and a Fourier transformation with respect to the 
spatial coordinates x A tangent to the boundary (2 ^ A ^ n), i.e., we write u as 
a superposition of modes 

u(t, x l ) = ^(x 1 ) exp(st + iujax a ) (43) 

with s E C, Re(s) > and u>a G K- By inserting Q43JI into the evolution equations 
H38f) . we obtain the Fourier-Laplace-transformed version 

su = A 1 d 1 u + iuj A A A u. (44) 

Partitioning the matrix 

sI m -iu, A A A = B=( I 12 ) (45) 

in the same fashion as in (|41(l (I m being the m x m unit matrix), equation l|44(l splits 
into 

=B n z + B 12 v, (46) 
A^i; = B 12 T z + B 22 v. (47) 

Because the A 1 are symmetric, the matrix B and in particular B\\ in l|45ll has 
only nonzero eigenvalues. Hence we may use the algebraic relations l|46() in order 
to eliminate z in favour of v. We are left with the system of ordinary differential 
equations (ODEs) 

dtv = Mv, M = A- 1 (B 22 -B 12 r B 11 - 1 B 12 ). (48) 
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A simple argument [341121) shows that M has precisely I eigenvalues Ai, A2, • ■ • , A; 
with negative real part. Let W\, w 2 , ■ ■ ■ , w ; denote the corresponding eigenvectors. 
Then the general L 2 solution of l|48|l is given by 

1 

v(s, a; 1 , = 2^ a j w j( s i u ) exp(A J -x 1 ) (49) 
i=i 

with complex integration constants o~j . (If the eigenvectors do not span the space then 
some of the aj have to be replaced with polynomials in x 1 - this will occur as a special 
case in our application.) 

The constants aj are determined by the boundary conditions: substituting <|49ll 
into the Fourier-Laplace transform of (|40|l . 

S 1 diu + \u A S A u = sh, (50) 

we obtain a system of linear equations 

C{a,w)<r = h, (51) 

where C is an / x I matrix. 

Consider first the homogeneous boundary conditions, h = 0. Suppose that 
detC = for some s with Re(s) > 0. Then l|51|l has a nontrivial solution, and 
hence the IBVP (|38H40[) with homogeneous boundary conditions has a non-trivial 
solution of the form l|43[) . Now observe that with 143[) . 

u(t, x l ) = ufax 1 ) exp(asi + \aujAX A ) (52) 

is also a solution, for any a > 0. Thus we can find solutions which grow exponentially 
at an arbitrarily fast rate, and the IBVP is ill-posed. Such solutions are called strong 
instabilities. (In a numerical simulation, a is determined by the highest frequency 
that can be represented on the grid. Hence the growth rate of the instability increases 
with resolution.) We conclude that the determinant condition 

det C{s,uj) 7^0 for Re(s) > (53) 

is a necessary condition for well-posedness. 

Next, we consider the inhomogeneous boundary conditions 1)40(1. Formally, we can 
solve l|51|l for the integration constants <r provided that the determinant condition is 
satisfied. What remains to be shown is that the solution l|49|) can be bounded in terms 
of the boundary data, 

\v(s,0,uj)\^K\h(s,uj)\, (54) 

with a constant K > that is independent of s and oj. This is known as the Kreiss 
condition [21] ■ Provided that the eigenvectors in (|49|l are normalized in such a way 
that they remain finite as Re(s) J. and as \s\ — > 00, l|54|) is equivalent to 

detC(s,w)^0 for Re(s) ^ 0. (55) 

Comparing this with (|53|l . the additional requirement is that there be no zeros s of 
det C with Rc(s) = 0. Such zeros are known as generalized eigenvalues. 

If the Kreiss condition is satisfied then it follows immediately that the IBVP l|38l 
I40|) is well-posed for vanishing initial data. The unique solution is given by the inverse 
of the Fourier-Laplace transform, which is well-defined because of the bound (|54|) . In 
particular, we obtain an estimate of the form 

[ \\u(;T)\\ 2 Q dr4K T [ \\h(;T)\\l Q dT (56) 
Jo Jo 
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in every time interval ^ t ^ T, where the constant Kj- is independent of the data h 
(the norms are L 2 norms over the half-space and the boundary, respectively) . We say 
that the system is boundary- stable |2U) . It is this definition of stability that we shall 
prove for the generalized harmonic evolution system (section 13. 2|) . 

Closely related to boundary stability is the concept of strong well-posedness in 
the generalized sense |21l 1221 1277) . This requires that if we add a source term F(t, x 1 ) 
to the right-hand-side of the evolution equations (|3"%)) , the following estimate holds, 




(recall that v refers to the modes with non-zero speeds). The crucial ingredient in 
proving (j57(l given Ij56|l is the construction of a symmetrizer |35U34ll2T)] . In doing this, 
it is usually assumed that the boundary conditions are in maximally dissipative form 



v-=Sv++g, (58) 

where v~ and v + are the ingoing and outgoing (non-zero speed) modes and S is a 
(for our purposes constant) matrix. The symmetrizer method docs not appear to be 
applicable if the boundary conditions are of the differential type (|4*n)l . 

So far we have assumed that the initial data / vanish. One can always treat the 
case of general initial data by considering e.g. u' = u — e - */, so that the problem for 
ul has zero initial data. However, in the evolution equations for it' there will be an 
additional source term containing derivatives of /, which will appear in the estimate 
generalizing (|57|l . What one would like instead is an estimate of the form 




without any derivatives of /. This is referred to as strong well-posedness. Majda 
and Osher [31] showed in the present case of a uniformly characteristic boundary 
that the Kreiss condition is also sufficient for strong well-posedness provided that 
the boundary conditions are maximally dissipative. However, nothing is known in 
general for boundary conditions of the form 1|40JI . In section 0] we will argue that the 
Kreiss condition is still useful in order to rule out certain weak instabilities that have 
a polynomial time dependence. 

We close this section with some remarks on the scope of the Fourier-Laplace 
method. Technically, one can only apply Fourier and Laplace transforms if the 
evolution equations and boundary conditions have constant coefficients. However, 
well-posedness is a concept that is associated with the behaviour of the high-frequency 
components of the solution. In many cases one can show that a problem with variable 
coefficients is strongly well-posed if all problems obtained by freezing the coefficients 
are strongly well-posed ([22 Theorem 8.4.9], see [2Fj| for the extension from strictly 
hyperbolic to symmetric hyperbolic systems, see also the discussion in [201 ) - Strong 
well-posedness can be further extended from linear to quasilinear systems (such as the 
Einstein equations), in which case the estimates will in general only hold in a finite 
time interval (see for example [23 sec. 8.5], [7| [S])- More general spatial domains 
Q, than the half-space can easily be treated by splitting <9f2 into a finite number of 
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portions, each of which can be smoothly mapped to the half-space sec. 9.6.2]. 
One should note, however, that the proofs of the above extensions usually assume the 
boundary conditions to be maximally dissipative - the situation is much less clear for 
differential boundary conditions. 



3. 2. Application to the GH system 

In this section, we apply the Fourier-Laplace technique to the generalized harmonic 
system of section |3 For previous applications of this method to various formulations 
of the Einstein equations, see e.g. |51 12H |2"1 H I231 |2"U| . 

First we construct the most general frozen-coefficient problem by considering 
the limit of high-frequency perturbations of an arbitrary fixed background spacetime. 
Consider a point p at the boundary. By rescaling and rotating the spatial coordinates 
if necessary, we can achieve that the 3-metric at p is gij — &y , and by rescaling the 
time coordinate, we may also assume that the lapse at p is N — 1. However, we have to 
allow for an arbitrary shift vector. (If we performed a coordinate transformation that 
affected the shift vector, we would have to consider a moving boundary.) Furthermore, 
we may assume that the boundary is located at x 1 = x — and that the domain of 
integration is x > 0. 

Hence the evolution equations become 

dtipab = N A d A 4> ab , (60) 

d t U ab = N x d x U ab - d k $ kab - l2 N k d k ^ ab , (61) 

d t $jab = N x d x $ jab - djli ab + j 2 djtpab, (62) 

where we have introduced the operator dt = dt — N a 3a, N l refers to the (constant) 
background shift vector, and it is now understood that spatial indices are raised and 
lowered with the unit metric. 

We consider an ansatz similar to (|43ll . 

u(t, x l ) = u(x) exp [st + ilo a {x a + N A t)] (63) 

for all the dependent variables u = {i/j ab , H^, $i a f>}, where s £ C with Re(s) > 
and x A = {y, z}. Because of the rotational invariance of the evolution equations and 
boundary conditions, we may without loss of generality take loa — uj5a v , say, with 
uj > (note that uj = can be excluded because we have taken the high-frequency 
limit). It is convenient to eliminate w from the following equations by defining £ = lux 
and C = s/ui. 

Let us first consider the case N x — 0. Inserting the ansatz 1|63J) into the evolution 
equations l|6UH62|l . we first obtain the algebraic relations 

C^Pab = iN y i> ab , (64) 

($ yab = -ifl ab , (65) 

(® Z ab = 0. (66) 

From these we deduce that ip ab — $ zab = 0, and we eliminate <& ya b m favour of H ab - 
The remaining two evolution equations are turned into the ODE system 

*(&)-(-/-* ?)(*-)■ 

The matrix has eigenvalues 



\f = ±v/l + C 2 (68) 
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(the branch of the square root is chosen such that Re(A^) > for Re(C) > 0). Since 
we are only interested in L 2 solutions, i.e., solutions that decay as x — > oo, we have 
to pick the eigenvalue A^ . The corresponding eigenvector is (£, — A^~) T and hence the 
general L 2 solution is 



( C \ 



0~lab e 



a; 



c 

-K 

— i 




(69) 



§xab 
§yab 

V J V J 

with arbitrary complex constants o~x a b- The normalization constant K\ (and all similar 
constants in the following) is equal to the norm of the vector immediately following 
it. 

Next we discuss the case = N x ^ 0. We obtain tp a i, = as before, but now the 
ODE system reads 

,,2/ 



0, 



( \ 

fxafe 
%ab 
\ ®zab J 

where we have introduced the shorthand j 2 



( 



-i 2 K 

-7 2 C 

i/0 





-7 2 C 






-i7 2 
t/0 








C/4 
(i 





( Uab \ 




$>xab 




%ab 


) 


\ ®zab J 



Af = 7 2 (-/3C± VC 2 +7" 2 ), 



-0 2 ) 
A 2 = 



The eigenvalues are 



c 

0' 



For > 0, only Aj has negative real part (recall that we are assuming 
and the corresponding solution is 

/ _ 



%ab 
\ $zab ) 



o- lab e"^ k i 1 



( C-0K \ 
-XT 



V 







(70) 

(71) 
<1) 

(72) 



/ 



One observes that this contains the above result Ij69|) for — in a regular way so 
that we do not need to discuss these two cases separately in the following. 

For < 0, the eigenvalue A2 also has negative real part and hence the general L 2 
solution is 

( C - 0K \ 

-XT 



( n afc 

$xab 
®yab 
\ $zab 



\ 



0~lab<S 



AT 



^r 1 



V 



3 A 2 « 





/ 



C lab Kn 



J 
\ 



V 



(73) 



+ rr 

£ +&3ab q 

/ V 1 J 

More care is needed if £ = — because in that case, A^ = A2 = — 1 and the 
corresponding eigenvectors are degenerate (this is the only case in which such a 
situation occurs). The most general solution is now of the form 



( n " b \ 






f -P \ 




= e-« 


0~lab K 1 1 





V $zab J 






V ) 



( 






( \ 




1 















+ 0~3ab 







V o 1 


J 




V 1 J 
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(74) 

In order to work out the Fourier-Laplace transform of the boundary conditions, 
it is useful to note that in our set-up, {t a ,(d/dx) a — —n a ,(d/dy) a 7 (d/dz) a } form 
an orthonormal basis, and Pij = S?Sj + 8*6*. An index t will be used to denote 
contraction with the timelike normal t a . We remark that the c'iV'ab terms hidden in 
the omitted Cj a & terms in <|14l 1151) do not enter in the following because we already 
know that ip ab = 0. 

The Fourier-Laplace transform of the constraint-preserving boundary conditions 
l|2(ill then reads 

ChP = c°- = d i (-|n s - n xi - \i\ xx - ±n ra - ±n Z2 

IKtt ~ Kxt ~ \Kxx - \Kyy - hKzz) (75) 



2 xtt ^xxt 2 xxx 2 X VV 2 

+ 1 ~ ~ \Kxx ~ \Kvy ~ \K»* - Kyi 

Chf ] = c° x - = d t (-|n s - u ix - ±n ra + ±n ra + \t\ zz 

_I$ $ . _ 1$ + I<t> 4- la ^ (761 

2 xti Y a:ta 2 111 T 2 T 2 Izz J ^ 1 "/ 

+ i (-flyx - $ yix - Kyx) , 



1 - 



2 it ^ 2 xx 2 W 2 zz C^) 



C^i C) = eg" = di ^il iz -fl xz -^ xiz -^ xxz )+i[~fl yz -^ yiz -^ xvz ). (78) 
In the f3 < case, we have the additional constraint-preserving boundary conditions 

CC* = 5 "^ 6 = dzKab ~ iKab, (79) 

(h ( 6 % = ~c xzab = d^ zab . (80) 

Although it has to vanish for a solution that satisfies the constraints, we have added 
boundary data to all the constraint-preserving boundary conditions in order to 
account for the inevitable numerical truncation error. 
The physical boundary conditions (|31|l become 

= Wyy = 8 S (-^ + ^-1$^ + ^^) (8l) 

+ i (iftyx + Kxy + Kzy ~ \®y XX ~ ^Kyt + \® yxt) . 

(h ( p=w- z = d^(-fl yz -^ xyz ) (82) 

±ri Z * + Kxz + \§zzz - \®zxx - \®zyy - \§ xzi + ^ zxi ) ■ 

Note that the terms involving a normal derivative of v? ab = ijj ab in 1-iL'l) do not 
contribute because the Fourier-Laplace transform of "006 vanishes as a consequence 
of the evolution equations. 
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Finally, the gauge boundary conditions l|3t)|) are 

±fL u -fl xt + $ xii -$ xxi , (83) 
hf ] = % " n ra + $ x£c - $ xxx , (84) 

>4 G) = % - n x , + <b xiy ~ <s> xxy , (85) 

= U iz - U xz + $ xiz - $ xxz . (86) 

The next step in the calculation is to insert the general L 2 solutions into the 
Fourier-Laplace transforms of the boundary conditions. We obtain a system of linear 
equations 

C(C)<r = h. (87) 

We need to evaluate the determinant of the matrix C . The resulting expressions are 
rather lengthy and some computer algebra is helpful here. A programme (available 
from the author upon request) has been written in the computer algebra language 
REDUCE [37J in order to obtain the following results. We discuss the cases (3 ^ 
and (3 < separately. 

Suppose first that (3 ^ 0. In this case the general solution is l|72ll and the boundary 
conditions are l|75H78l l5TH8l)}l . We find 

^±™r m 

The function £ — > £ — (1 + (3)X[ maps the right half of the complex plane to the right 
half of the complex plane minus a circle of radius [(l + /3)/(l — /3)] 1 / 2 . Hence detC has 
no zeros £ with Re(£) > and is bounded away from zero even as Re(£) — > 0. In the 
limit |C| — > oo, we have |A^| ~ |£|, \K X \ ~ |£| and hence detC is bounded away from 
zero, too. We conclude that both the determinant condition and the Kreiss condition 
are satisfied. 

Next, we consider the case (3 < 0. Suppose first that £ ^ —(3 so that the general 
decaying solution is (|73|l . Plugging this into the boundary conditions l|75H8fijl yields 

-[C- (l+y9)Ar] 16 (C 2 ~ (3 2 ) w 
8( 16 f3 20 Kl K? 

By the same argument as above, detC has no zeros £ ~P with Re(£) ^ 0. In the 
degenerate case £ = —(3, we have to use the solution Q74J1 instead and find 

(K A- K \ w 

We conclude that for arbitrary shift at the boundary (with subluminal normal 
component), our initial-boundary value problem satisfies the determinant condition 
and the Kreiss condition in the high-frequency limit. In particular, there are no strong 
instabilities. 

It is easy to see that this result is unchanged if the physical boundary conditions 
13111 are replaced with the improved conditions (|34|l of [33]. In the frozen-coefficient 
approximation used in this section, they correspond to successive applications of the 
operator (t a + n a )d a = dt — (1 + (3)d x to the boundary conditions J3jJ. Each produces 
an additional factor of £— (1 +j3)X[ in the expressions for det C, and as argued above, 
this factor is bounded away from zero. 



detC " 0^16^20 ylOylO ■ 
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Finally, we ask what happens if we replace the gauge boundary conditions (|35[1 
with the alternative conditions (|37[1 . In this case, we obtain 

■ -[c-(i + /?)Ar] 8 (c-/3Ar) 8 



det C = I 



-[C-(i 



8( 

-P)K 



if P > 0, 



if p < 0, 



(91) 



8(16/320^1 

so that det C always has a zero 

C = ~P- (92) 

Both the determinant condition and the Kreiss condition are satisfied if the shift 
points towards the interior at the boundary (p > 0). However if the shift points 
towards the exterior (P < 0), the IBVP is ill-posed. For tangential shift (P = 0), the 
determinant condition is satisfied (and hence there are no strong instabilities) but the 
Kreiss condition is violated (C = is a generalized eigenvalue). What happens in this 
case will be investigated in the following section. 



4. Weak instabilities 



As we have seen in the previous section, the determinant condition i|53|) is sufficient 
to detect strong instabilities of the form l|43l) . The question remains whether a given 
IBVP admits ill-posed modes with "milder" than exponential time dependence, so- 
called weak instabilities. Examples of such instabilities with linear time dependence 
have been found in a formulation of Maxwell's equations in a gauge with vanishing 
electrostatic potential |2H and in an Einstein-Christoffel formulation of general 
relativity 0]. 

In this section, we develop a technique that can be used to systematically search 
for weak instabilities with polynomial time dependence. The criterion turns out to be 
closely related to the Kreiss condition. We then apply our method to the generalized 
harmonic system and show that it does not suffer from such instabilities. Furthermore, 
we show that if we use the alternative gauge boundary conditions (|37|1 instead, a weak 
instability does exist, and we construct it explicitly. 



4-1- General theory 

We consider again the general IBVP 138H40fl . As before, we Fourier-transform with 
respect to the coordinates x A tangential to the boundary. However, instead of 
assuming an exponential time dependence as in the Laplace transform (|43|) . we now 
look for modes whose time dependence is given by a polynomial of order p ^ 1, i.e., 
we make the ansatz 

u{t,x l ) = e-Kv(iLuCj A x A )Y uMfw 1 )^, (93) 

where to > 0, u>a& a = 1, and ^ 0. Substituting this into the evolution equations 
H38fl and evaluating order by order in t, we find (setting £ = lux 1 ) 

A l d^u {p) + \& A A A vP> = 0, (94) 

A l d i u^ p - l) + ilu a A a u < -p-^ = fiW, (95) 

A 1 d l .u (u) +iu) A A A u {u) =u {,J+l) , 0^v^p-2. (96) 
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Let us also expand the boundary data h in a similar fashion as in (|93|) . up to 
order p — 1 for a reason that will become clear shortly, 

fc(i,:r A ) = exp(iLj(Jj A x A )Y hM^¥-. (97) 

i/=0 

The boundary conditions 140(1 then read 

S^fiW + ito A S A u {p) = 0, (98) 

S 1 3 c u< p - 1 ) + iwx^ii^- 15 = 0, (99) 

fi^fiM + iw A 5 A tt (l/) = h(" +1 > , < i/ s$ p - 2. (100) 

Suppose that a solution tt = (u^ ,u^\. .., ) satisfying equations (|94II96() and 
l|98H100|l exists. Then the solution at times t > is of order 0(ui p ) whereas the initial 
data is O(w ) and (the time integral of) the boundary data is 0(w p_1 ). Hence the 
solution cannot be estimated in terms of the initial and boundary data as ui — > oo and 
the IBVP is ill-posed. We may actually argue that it is already ill-posed if a (nonzero) 
solution of (JMUnil) exists that obeys the top two boundary conditions (|9"5H9"9"|) , for we 
can always choose the boundary data h such that the remaining boundary conditions 
H 100(1 are satisfied. 

Let us look more closely at equations (19411 and ((98(1 : these are identical with 
equations l(44|) and ((5011 in section ETTl for s = 0. We conclude that if a weak instability 
exists, s — is a generalized eigenvalue. (Strictly speaking, one has to consider the 
Fourier-Laplace solutions for Re(s) > and take the limit s — > 0. We assume here 
that the solutions are continuous at s = 0, which is the case in all the examples we 
discuss.) Therefore, if the Kreiss condition is satisfied, there are no weak instabilities 
of any polynomial order p. 

On the other hand, s = being a generalized eigenvalue does not automatically 
imply that there is a weak instability. One still has to solve the coupled ODE system 
l(94H96(l and impose the boundary conditions l(98H99(l on its solution. 

In [23 it is claimed that (for a certain choice of parameters) the Maxwell system 
considered there admits a weak instability (with p = 1) even though the Kreiss 
condition is satisfied. This is obviously a contradiction to the result that we have 
just derived. However, the eigenvectors Wj appearing in the general L 2 solution in 
(241 eq.(51)] are normalized in such a way that they become infinite as s — * 0, and 
the estimate (154(1 does not obtain. If the correct normalization is used, one finds that 
s = is indeed a generalized eigenvalue. 

4-2. Application to the GH system 

To see how the method outlined above works in practice, we apply it to the generalized 
harmonic Einstein equations in the high-frequency limit (equations l(60H62(l V 

First, we take the gauge boundary conditions to be 1(85(1 . We have already seen 
in section 13.21 that the Kreiss condition is satisfied, and hence there are no weak 
instabilities with polynomial time dependence. (One can also carry out the explicit 
calculation below, obtaining the same result.) 

Next, we consider the alternative gauge boundary conditions ((37(1 . and we take 
the shift to be tangential at the boundary (j3 — 0). As discussed in section ET2l the 
Kreiss condition is violated in this case, with s = being a generalized eigenvalue. 
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This indicates that there might be a weak instability. Let us search for one with p = 1 
and homogeneous boundary data. Equation (|94|) implies 

-9^-1^ = 0, (101) 
-ifi£ = 0, (102) 



and H95|) reads 



-^=*a (104) 

-in^ = (105) 

= CI (106) 



Substituting l|lU4llU5f> into (|TUl"|) . we obtain 



^fiS^nS =► ni° 6 3 =^e-«. (107) 

Assuming that Ci a b = initially, i.e., 

*s=wff, €i=o, (los) 

equation (|103fl together with (|102fl implies similarly 

<W=$? =► tfff = «-W«. (109) 
Hence the most general L 2 solution (consistent with Ci a b = 0) is of the form 

(ni° b \ CI = e "' (^o6. -^6, i<T2a6, 0), (110) 

(niM'i CV) = e"«(0, <7 lab , -i<r lo6) 0). (Ill) 

Only the mode associated with o~i a b corresponds to a solution with ^ 0. 

Substituting this mode into the Fourier transform of the homogeneous boundary 
conditions and evaluating order by order in t, one obtains a 20 x 10 linear system 
of equations for the constants o~i a b- The only nontrivial solution is given by 

o-ixx = 1, vxxy = — i, o- lyy = — 1, all other a kah = 0. (112) 

It corresponds to setting 

riy = f,ij, n 0a = 0, 

®kij = ~tf,kij, $fc0a = 0, (113) 

which is of the same type as the weak instabilities found in [2U E] . 

A similar calculation shows that there are no weak instabilities of order p 2 for 
homogeneous boundary data. However, as explained in section T4. II this is not true for 
inhomogeneous boundary data. 



Stable radiation- controlling boundary conditions 



18 



5. Robust stability tests 

In this section, we perform a series of robust stability tests on our numerical 
implementation of the GH system. This test involves adding small random 
perturbations to a known exact solution. Thus all possible frequencies at a given 
resolution are excited. By increasing the resolution and looking at the growth rate of 
the numerical solution, this test is very effective in spotting ill-posed modes. The 
robust stability test was suggested in the context of the Cauchy problem by the 
"Apples with Apples" collaboration see also for an application to the spectral 
evolution code used for this work). Previous studies in which a version of the test was 
used in relation to stability of the IB VP include |3U1 IT71 IT%1 H HI ) . 

For most of the tests, the background solution is taken to be flat space with a 
constant shift N l , 

ds 2 = -dt 2 + Sij {dx l + N i dt)(dx j + N j dt). (114) 

We remark that this is the most general background solution in the high-frequency 
limit, as discussed in section EOl The spatial domain is taken to be of topology T 2 x R. 
This domain is intended to be as simple as possible, avoiding additional complications 
due to sharp corners and edges. We take x,y,z G [—0.5, 0.5]. The boundary conditions 
discussed in this paper are imposed at x = ±0.5, and periodic boundary conditions 
are imposed in the y and z directions. 

For the last test, we consider instead the Schwarzschild solution in Kerr-Schild 
coordinates, 

ds 2 = - ( 1 - —) dt 2 + —dr dt+(l+ —) dr 2 + r 2 (d8 2 + sin 2 9d<j> 2 ) . (115) 



r J r 

In this case, the spatial domain is taken to be a spherical shell extending from 
^min = 1-8M (just inside the event horizon) to r max = 11. 8M. The boundary 
conditions discussed in this paper are imposed at the outer boundary r = r max . No 
boundary conditions are imposed at the inner boundary r — r min . See also for 
a spherically symmetric robust stability test on Schwarzschild spacetime in Painlevc- 
Gullstrand coordinates. 

As initial data, we compute the variables tp a b, n a & and <&i a b corresponding to the 
exact solutions ( jl 141 1115)1 and add random noise of amplitude e = 10 -10 (in the flat- 
space case) or e = 10 -6 (in the Schwarzschild case) to them. In addition, random noise 
of the same amplitude is added both to the right-hand-sides of the expressions for the 
time derivatives at the boundary (|118fl and to the right-hand-sides of the evolution 
equations. In this way, we probe the influence of nonzero initial and source data as 
well, which could not be analyzed using the methods of section 13) 

For more physically realistic tests involving black hole spacetimes, we refer the 
reader to |27j . In particular, an ingoing gravitational wave was injected through the 
outer boundary (i.e., data M p ) were prescribed in equation Q31|)'l. The gravitational 
radiation scattered back by the black hole was extracted, correctly reproducing the 
dominant quasi-normal mode oscillation. 



5.1. Numerical implementation 

We use pseudospectral methods as described for example in 

For the flat-space tests on T 2 x R, the numerical solution is expanded into 
Chebyshev polynomials in x and into Fourier series in y and z. We use between 
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9 and 27 basis functions in the x and y directions (these are typical resolutions in 
realistic numerical relativity simulations using spectral methods) and 3 basis functions 
in the z direction. Reducing the effective dimension by one enables us to evolve 
up to final times of t = 1000 within a tolerable runtime of the code even for the 
higher resolutions. We have checked for moderate resolutions that the behaviour is 
qualitatively unchanged in fully three-dimensional runs. A fourth-order Runge-Kutta 
scheme is used for the time integration, with a Courant factor of Af/Ax m in = 1-5, 
where Xmin is the minimum spacing between the (in the Chebyshev case, non-uniform) 
pseudospectral collocation points. 

For the Schwarzschild test, the numerical solution is expanded into Chebyshev 
polynomials in r and spherical harmonics in the angular directions. Here we use 
between 15 and 27 radial and between 11 and 19 angular basis functions. The same 
time integration scheme with the same Courant factor is used as in the flat-space case. 
The highest-order spherical harmonics are filtered as described for example in |3] . No 
filtering is applied to the Chebyshev and Fourier expansions. 

The boundary conditions are implemented by prescribing the time derivatives of 
the characteristic variables at the boundary 42 . For the Neumann- like constraint- 
preserving boundary conditions ((26129(1 and physical boundary conditions 1|> . we 
use the evolution equations for the incoming characteristic variables in order to treat 
normal derivatives for time derivatives, 

d t u ab = (N n + N)d n u 1 ab + tangential derivatives, (H6) 

dtu\ ab = N n d n u 2 Aab + tangential derivatives. (H7) 

(The notation d instead of d indicates that the partial derivative acts only on the 
fundamental variables, not on the eigenvectors needed to form the characteristic 
variables |27|.') For the gauge boundary conditions, we simply take a time derivative 
of 135(1 . In this way, we obtain expressions for d t u 2 , P^dtu 1 -, pWdtu 1 -, and 
P^dtu 1 ^ (omitting the spacetime indices for simplicity). These are then combined 
to form 

dtu 1 - = P^dtu 1 - + P^dtu 1 - + P^dku 1 - (118) 

at the boundary (note that p( c ) , p( p ) and p( G ) are mutually orthogonal projection 
operators adding up to the identity). In order to implement the alternative gauge 
boundary conditions (|37() . we first fill P^dtu 1 - and P^d t u 1 as before and then 
set 

P^dtu 1 - = QP^dtu 1 - + (Q- P^)d t (u 1+ + 2 72 V>), (119) 

where the operator Q has the properties P^Q = P^Q = 0, P^Q = Q, 
QP(C) = Q, Qp( p ) = QP( G ) = 0, and t ■ (P^ + Q) = 0. It then follows from 
H118|l and 1(1 0(1 that dtit • II) = at the boundary, as desired. An explicit expression 
for Q is given by 

Qab cd = k a k b k c k d + 2k [a P b) ( c k^ - 2k (a l b) k c k d . (120) 
5.2. Numerical results 

The following plots show the L 2 norms of the error 

£ = ViWab) 2 + (Sn ab )2 + (5$ ia6 )2 (121) 

and of the constraints 

C = y/(C a ) 2 + (T a Y + (C w ) 2 + (C lab f + (C ijab ) 2 (122) 
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Figure 1. Robust stability test on the flat-space background 11141 with 
constraint-damping. L 2 -norms of the error (left) and the constraints (right) for 
two different shift vectors N 1 and for four different resolutions M = Afx = My ■ 



as functions of time. In the above, the notation (M a fe... c ) 2 refers to the sum of the 
squares of the Cartesian components of M a b---c- The differences S in (|1 2 1|> are taken 
with respect to the background solutions (jl 141 1115(1 . When computing the quantity 
Sipab, we subtract its spatial average from it. This is because the constant mode of 
ipab typically has a large linear drift in time caused by a non-zero constant mode of 
n a fe (which is not eliminated by the differential boundary conditions). Note that this 
procedure does not affect the higher-frequency modes. 

We begin with the flat-space background (|114fl . We run the test for two different 
shift vectors, N l — (0,0,0) and N' 1 = (0.5,0.5,0). In the latter case, the shift points 
towards the interior at the x — —0.5 boundary (corresponding to the (3 > case in 
section and towards the exterior at the x = 0.5 boundary (corresponding to the 
(3 < case), and the shift also has a component tangential to the boundary. Figure 
^ shows that in both cases, the error and the constraints remain of the same order as 
initially up to t = 1000 (and presumably forever). For this run we included constraint 
damping with parameters 70 = 72 = 1 (the same choice was made in |27| in order to 
obtain long-term stable black hole evolutions) . This leads to the sharp initial decrease 
of the constraint violations (note that the randomly perturbed initial data do not 
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Figure 2. Robust stability test on the flat-space background 11141 without 
constraint damping. L 2 -norms of the error (left) and the constraints (right) for 
two different shift vectors N 1 and for four different resolutions M = Afx = My ■ 



satisfy the constraints). Because one might suspect that the constraint damping might 
somehow hide potential instabilities, we rerun the test without constraint damping 
(figure EJ • Now the constraints grow slightly (as expected) but still there is no sign of 
an instability. 

Next, we turn to the Schwarzschild background 1)115(1 . Here we choose the 
amplitude of the random perturbations somewhat higher than in the flat-space case 
(10 -6 rather than 10~ 10 ) so that it is much larger than the error incurred during an 
evolution of the unperturbed Schwarzschild spacetime for the resolutions considered. 
In this case, the constraint damping is essential in order to avoid exponential growth 
of the constraints 27 . Figure [3] shows the results of this test, using both the I? 
norm and the L°° norm (which is more sensitive to local effects on the boundary). 
There are now significant oscillations in the error due to the much higher amplitude of 
the perturbations but they grow only linearly on average and at a rate that does not 
appear to increase significantly with resolution. The constraints remain essentially 
constant thanks to the constraint damping. We remark that as a consequence of 
the stability analysis in section T6. 21 flat space and Schwarzschild spacetime should be 
equivalent with regard to stability in the high frequency limit. 
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The numerical results strongly suggests that the initial-boundary value problem 
is indeed well-posed, even if nontrivial initial and source data are included. Previous 
studies [41)1 1171 1181 0J presented similar numerical evidence (for different 
formulations), although only initial and boundary data but not source data were 
included in those robust stability tests. 

Finally, we return to flat space i|114l) but now we replace the gauge boundary 
conditions l|35() with the alternative set (|37|l . The top half of figure 0] shows the 
results of the robust stability test for shift vector N l = (0.5,0,0). Now the error 
grows exponentially at a rate that increases with resolution (note also the timescale as 
compared with hguresmilj!. This is what we expect because the initial-boundary value 
problem is ill-posed in this case (section [3. 2fl . The bottom half of figure 21 shows what 
happens if we vary the normal component N x of the shift. As predicted by the Fourier- 
Laplace analysis in section l3~21 fcf. equation H92fl ). the growth rate depends linearly 
on N x . Numerically, we find a nonzero offset of the growth rate at N x — 0, which 
is not obvious from the analysis. We remark however that because of the random 
perturbations, N x will never exactly vanish at the boundary, and furthermore the 
argument in section0]indicates that even for N x = 0, there are polynomial instabilities 
of arbitrarily high polynomial order if nonzero boundary data are included. 
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Figure 4. Robust stability test on the flat-space background 11141 , using the 
alternative gauge boundary conditions 1371 , Top half: L 2 norm of the error 
for four different resolutions J\f = Afx = Afy, using N x = 0.5 for the normal 
component of the shift. Bottom half: L 2 norm of the error for various values of 
N x at fixed resolution J\f = 9, and linear fit to the growth rate. (In all cases, 
N« = N z = and 70 = 72 = 1.) 



6. Conclusions 

We considered the initial-boundary value problem for a first-order formulation of the 
Einstein equations in generalized harmonic gauge. This system was derived in |27| 
and has proven very successful in obtaining long-term stable black hole evolutions. 
The boundary conditions we considered have the special property that they control 
the incoming gravitational radiation via the incoming fields of the Weyl tensor 
0121 El El- We believe that this is essential in order to obtain reliable information 
about the gravitational radiation emitted from a compact source if the domain of 
integration has artificial timelike boundaries, as is a common situation in numerical 
relativity. In addition, the boundary conditions eliminate the incoming constraint 
fields [51 1101 ITU 13*1 114 ) . in which we believe they are superior to constraint-preserving 
boundary conditions of Dirichlct type |2l)| from a numerical point of view. 

In section"""] we analyzed the well-posedness of the initial-boundary value problem 
using the Fourier-Laplace technique """J [313 H21 [ffl HH1 1201' This required taking 



Stable radiation- controlling boundary conditions 



24 



the high-frequency, or frozen-coefficient approximation. To allow for an arbitrary 
background spacetime, we had to take into a account an arbitrary shift vector at 
the boundary. This generalizes the result stated in |27], where only a tangential 
shift was considered. We showed that the Kreiss condition is satisfied, which implies 
that the system is boundary-stable (i.e., the solution can be estimated in terms of 
the boundary data). Unlike for maximally dissipative boundary conditions, it is not 
known in the present case of differential boundary conditions whether (or under which 
additional assumptions) the Kreiss condition is also sufficient for well-posedness if non- 
trivial initial and source data are included. It would be of considerable interest to the 
numerical relativity community to obtain a general theorem covering this case. 

It has been claimed [^] 0] that systems with differential boundary conditions 
might admit weak instabilities with milder than exponential time dependence even 
if the Kreiss condition is satisfied. In section we considered instabilities with 
polynomial time dependence for a general first-order initial-boundary value problem. 
It was found that such instabilities are ruled out by the Kreiss condition. (More 
precisely, the condition is that be not a generalized eigenvalue.) For the generalized 
harmonic Einstein equations, it turned out that the choice of boundary conditions for 
the gauge degrees of freedom can be crucial for stability: for an innocent-looking set 
of alternative gauge boundary conditions that did not satisfy the Kreiss condition, 
we found a weak instability (with linear time dependence) if the shift was tangential 
at the boundary. However, as soon as the shift pointed towards the exterior, the 
weak instability was turned into a strong one with exponential time dependence - 
this demonstrates that taking into account a normal component of the shift can be 
important. 

Finally, we performed a numerical robust stability test |3B1 EH1 E3 E3 El El El 

14 1| . The background spacetime was taken to be either Minkowski space with a shift 
or Schwarzschild. We added small random perturbations to both the initial data, 
the boundary conditions and the right-hand-side of the evolution equations. The 
generalized harmonic evolution system performed very well on these tests, with the 
error and the constraints remaining of the same order of magnitude over 1000 light 
crossing times in the flat-space case, or 1000M in the Schwarzschild case. This strongly 
suggests that the initial-boundary value problem is likely to be well-posed even if 
nontrivial initial and source data are included. We also ran the test on the alternative 
set of gauge boundary conditions, finding resolution-dependent exponential growth 
of the error as predicted by the Fourier-Laplace analysis. The linear dependence of 
the growth rate on the normal component of the shift was also reproduced, although 
exponential growth was observed even in the limiting case of a tangential shift. 
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Appendix 

In this appendix, we prove a detail that we postponed in section l2~2l Recall that in the 
case of an outward-pointing shift ( N x < 0) , we could only impose constraint-preserving 
boundary conditions on the incoming constraint fields c° a ~ and c^ Aab . Here we show 
that this implies that the remaining incoming constraint fields c\ a and c\ Bab also 
vanish at the boundary, so that we have a full set of maximally dissipative boundary 
conditions for the constraint evolution system. We restrict ourselves to the high- 
frequency limit and use the Fourier-Laplace framework of section 13. 21 

Suppose that we impose the constraint-preserving boundary conditions proposed 
in section l2~21 in the N x < case: in Fourier-Laplace language, equations l|75H8t)|l with 
zero data h c . These imply a linear system 

C«<r = (A.l) 

for the constants o~i a b parametrizing the general Fourier-Laplace-transformed solutions 
(J73J). The kernel of C*W spans all solutions that satisfy our constraint-preserving 
boundary conditions. 

Next, we consider the remaining incoming constraints c\ a and c\ Bab . Their 
Fourier-Laplace transforms are given by 

& = di% xi + i(% yi + 1% + ifi xx + ±n TO + |n„), (A.2) 
<:.. =Sc(i*»»« + 5* w a-3*yi W -i*»«) + i(*i W x + n fe ), (A.3) 

<•:... =d^y X y+i(^yyy + ^ vU -^ yXX -^y ZZ +fl iv ), (A.4) 



u y t ~ w t,^yxt ' L \^yyt ' 2 xx tt ~ 2 XXxx ~ 2 XX VV ~ 2 

CyX ~ 2 •• ~I' V :>>' "~ 2" V U<<<! ~ 2 

^yy -"^ji!)TH2 t ji W t 71,,;, o -J.-,,,,-,, 2 ; ' '/-/ 

C 2 y Z = 9^ yxz + i{$ yyz + fifo), (A.5) 

c 2 

z 

7-2 



?i =d^ zxi + ^ zyh (A.6) 
cL = %{h*zxx + \® zi i - \®zyy ~ + M^x, (A.7) 

Cly = d^ zxy + i{\$ z yy + \* & ~ \§ zxx ~ (A.8) 

c 2 zz = d^ zxz +i^> zyz , (A.9) 

Cyzab = (A.10) 

Equations l)A.2HA.10|) imply another linear system for the integration constants ai a b, 

C {2) ct = Q. (A.ll) 

We need to check that any solution satisfying the boundary conditions (|75H86[1 
also satisfies equations <|A.2HA.10|) . in other words, that 

kerC (1) C kerC (2) . (A. 12) 

This is a purely algebraic condition, which is straightforward to check with our 
computer algebra programme. It is indeed satisfied. 
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